# clear working environment
rm(list = ls())

#### load the required packages ####
library(tidyverse)
library(PanelMatch)
library(countrycode)

#### load the clean data ####
df <- read_rds("data/dem-transitions-replication-data.rds") %>%
  glimpse()

#### make df into a properly formatted dataframe for PanelMatch ####

# begin coercing df into m
# df needs to be a dataframe class
m <- as.data.frame(df)
# check to see if m is a dataframe
class(m)

# country needs to be an integer to use PanelMatch
m$country <- as.integer(m$country)

# year needs to be an integer to use PanelMatch
m$year <- as.integer(m$year)

# create treatment variable
m$tr <- as.numeric(ifelse(m$v2x_regime >= 2, 1, 0))

# all other variables, except categorical variables, need to be numeric
m$mid <- as.numeric(m$mid)
m$fatal <- as.numeric(m$fatal)
m$ln_cinc <- as.numeric(m$ln_cinc)
m$ln_pecpp <- as.numeric(m$ln_pecpp)
m$terr_dispute <- as.numeric(m$terr_dispute)

# make sure all of the rows are unique
m <- distinct(m, country, year, .keep_all = TRUE)

# create "era" indicators
m <- m %>%
  mutate(era = case_when(year < 1946 ~ "Pre-Cold War",
                         year >= 1946 & year < 1990 ~ "Cold War",
                         year >= 1990 ~ "Post-Cold War")) %>%
  mutate(post_1945 = ifelse(year > 1945, "Post-1945", "Pre-1945")) %>%
  mutate(post_cw = ifelse(year >= 1990, "Post-Cold War", "Pre-/Cold War")) %>%
  glimpse()

# set seed to reproduce bootstrapped estimates
set.seed(1)

#### initialize a data frame for plotting the pooled outcome over three years ####
gg_three <- data.frame(outcome = c("All MIDs", "All MIDs", "Fatal MIDs", "Fatal MIDs"), number_periods = "3") %>%
  glimpse()

#### effect of democratic transition on MID participation (pooled three years) ####
# create panel data
panel <- PanelData(panel.data = m,
                   unit.id = "country",
                   time.id = "year",
                   treatment = "tr",
                   outcome = "mid")

# create matched sets
PMresults.mid.att <- PanelMatch(panel.data = panel, lag = 5, refinement.method = "CBPS.weight",
                                match.missing = TRUE, exact.match.variables = c("un_region"),
                                covs.formula = ~ I(lag(mid, 1:5)) + I(lag(ln_cinc, 1:5)) + I(lag(ln_pecpp, 1:5)) + I(lag(terr_dispute, 1:5)),
                                qoi = "att", lead = 1:3, forbid.treatment.reversal = FALSE,
                                use.diagonal.variance.matrix = TRUE)

# estimate the ATT
PEresults.mid.att <- PanelEstimate(panel.data = panel, sets = PMresults.mid.att, 
                                   se.method = "bootstrap", number.iterations = 10000,
                                   confidence.level = .90, moderator = "post_cw",
                                   pooled = TRUE)

#### effect of democratic transition on fatal MID participation (pooled three years) ####
# create panel data
panel <- PanelData(panel.data = m,
                   unit.id = "country",
                   time.id = "year",
                   treatment = "tr",
                   outcome = "fatal")

# create matched sets
PMresults.fatal.att <- PanelMatch(panel.data = panel, lag = 5, refinement.method = "CBPS.weight",
                                  match.missing = TRUE, exact.match.variables = c("un_region"),
                                  covs.formula = ~ I(lag(fatal, 1:5)) + I(lag(ln_cinc, 1:5)) + I(lag(ln_pecpp, 1:5)) + I(lag(terr_dispute, 1:5)),
                                  qoi = "att", lead = 1:3, forbid.treatment.reversal = FALSE,
                                  use.diagonal.variance.matrix = TRUE)

# estimate the ATT
PEresults.fatal.att <- PanelEstimate(panel.data = panel, sets = PMresults.fatal.att, 
                                     se.method = "bootstrap", number.iterations = 10000,
                                     confidence.level = .90, moderator = "post_cw",
                                     pooled = TRUE)

#### add estimates, and bootstrapped standard errors and confidence intervals to plotting data frame ####
gg_three$estimates <- c(PEresults.mid.att[["Pre-/Cold War"]][["estimate"]], PEresults.mid.att[["Post-Cold War"]][["estimate"]], PEresults.fatal.att[["Pre-/Cold War"]][["estimate"]], PEresults.fatal.att[["Post-Cold War"]][["estimate"]])
gg_three$boot_se <- c(PEresults.mid.att[["Pre-/Cold War"]][["standard.error"]], PEresults.mid.att[["Post-Cold War"]][["standard.error"]], PEresults.fatal.att[["Pre-/Cold War"]][["standard.error"]], PEresults.fatal.att[["Post-Cold War"]][["standard.error"]])
gg_three$boot_lo <- c(quantile(PEresults.mid.att[["Pre-/Cold War"]][["bootstrapped.estimates"]], probs = 0.05), quantile(PEresults.mid.att[["Post-Cold War"]][["bootstrapped.estimates"]], probs = 0.05), quantile(PEresults.fatal.att[["Pre-/Cold War"]][["bootstrapped.estimates"]], probs = 0.05), quantile(PEresults.fatal.att[["Post-Cold War"]][["bootstrapped.estimates"]], probs = 0.05))
gg_three$boot_hi <- c(quantile(PEresults.mid.att[["Pre-/Cold War"]][["bootstrapped.estimates"]], probs = 0.95), quantile(PEresults.mid.att[["Post-Cold War"]][["bootstrapped.estimates"]], probs = 0.95), quantile(PEresults.fatal.att[["Pre-/Cold War"]][["bootstrapped.estimates"]], probs = 0.95), quantile(PEresults.fatal.att[["Post-Cold War"]][["bootstrapped.estimates"]], probs = 0.95))

#### initialize a data frame for plotting the pooled outcomes over five years ####
gg_five <- data.frame(outcome = c("All MIDs", "All MIDs", "Fatal MIDs", "Fatal MIDs"), number_periods = "5") %>%
  glimpse()

#### effect of democratic transition on all MID participation (pooled five years) ####
# create panel data
panel <- PanelData(panel.data = m,
                   unit.id = "country",
                   time.id = "year",
                   treatment = "tr",
                   outcome = "mid")

# create matched sets
PMresults.mid.att <- PanelMatch(panel.data = panel, lag = 5, refinement.method = "CBPS.weight",
                                match.missing = TRUE, exact.match.variables = c("un_region"),
                                covs.formula = ~ I(lag(mid, 1:5)) + I(lag(ln_cinc, 1:5)) + I(lag(ln_pecpp, 1:5)) + I(lag(terr_dispute, 1:5)),
                                qoi = "att", lead = 1:5, forbid.treatment.reversal = FALSE,
                                use.diagonal.variance.matrix = TRUE)


# estimate the ATT
PEresults.mid.att <- PanelEstimate(panel.data = panel,
                                   sets = PMresults.mid.att, 
                                   se.method = "bootstrap", number.iterations = 10000,
                                   confidence.level = .90, moderator = "post_cw",
                                   pooled = TRUE)

# democratic transition on fatal MID participation (pooled five years)
# create panel data
panel <- PanelData(panel.data = m,
                   unit.id = "country",
                   time.id = "year",
                   treatment = "tr",
                   outcome = "fatal")

# create matched sets
PMresults.fatal.att <- PanelMatch(panel.data = panel, lag = 5, refinement.method = "CBPS.weight",
                                  match.missing = TRUE, exact.match.variables = c("un_region"),
                                  covs.formula = ~ I(lag(fatal, 1:5)) + I(lag(ln_cinc, 1:5)) + I(lag(ln_pecpp, 1:5)) + I(lag(terr_dispute, 1:5)),
                                  qoi = "att", lead = 1:5, forbid.treatment.reversal = FALSE,
                                  use.diagonal.variance.matrix = TRUE)

# estimate the ATT
PEresults.fatal.att <- PanelEstimate(panel.data = panel, 
                                     sets = PMresults.fatal.att,
                                     se.method = "bootstrap", number.iterations = 10000,
                                     confidence.level = .90, moderator = "post_cw",
                                     pooled = TRUE)

#### add estimates, and bootstrapped standard errors and confidence intervals to plotting data frame ####
gg_five$estimates <- c(PEresults.mid.att[["Pre-/Cold War"]][["estimate"]], PEresults.mid.att[["Post-Cold War"]][["estimate"]], PEresults.fatal.att[["Pre-/Cold War"]][["estimate"]], PEresults.fatal.att[["Post-Cold War"]][["estimate"]])
gg_five$boot_se <- c(PEresults.mid.att[["Pre-/Cold War"]][["standard.error"]], PEresults.mid.att[["Post-Cold War"]][["standard.error"]], PEresults.fatal.att[["Pre-/Cold War"]][["standard.error"]], PEresults.fatal.att[["Post-Cold War"]][["standard.error"]])
gg_five$boot_lo <- c(quantile(PEresults.mid.att[["Pre-/Cold War"]][["bootstrapped.estimates"]], probs = 0.05), quantile(PEresults.mid.att[["Post-Cold War"]][["bootstrapped.estimates"]], probs = 0.05), quantile(PEresults.fatal.att[["Pre-/Cold War"]][["bootstrapped.estimates"]], probs = 0.05), quantile(PEresults.fatal.att[["Post-Cold War"]][["bootstrapped.estimates"]], probs = 0.05))
gg_five$boot_hi <- c(quantile(PEresults.mid.att[["Pre-/Cold War"]][["bootstrapped.estimates"]], probs = 0.95), quantile(PEresults.mid.att[["Post-Cold War"]][["bootstrapped.estimates"]], probs = 0.95), quantile(PEresults.fatal.att[["Pre-/Cold War"]][["bootstrapped.estimates"]], probs = 0.95), quantile(PEresults.fatal.att[["Post-Cold War"]][["bootstrapped.estimates"]], probs = 0.95))

#### create plot of estimates ####
# combine into one data frame for plotting
gg <- rbind(gg_three, gg_five) %>%
  mutate(post_cw = c(rep(c("Pre-/Cold War", "Post-Cold War"), 4))) %>%
  glimpse()

gg$post_cw <- factor(gg$post_cw, levels = c("Pre-/Cold War", "Post-Cold War"))

# create Figure 9 in the main text
ggplot(data = gg, aes(x = outcome)) +
  geom_hline(aes(yintercept = 0), linetype = "dashed", linewidth = .25) +
  geom_pointrange(aes(y = estimates, ymin = boot_lo, ymax = boot_hi, shape = number_periods, color = number_periods), size = .4, linewidth = .70, position = position_dodge(width = .2)) +
  scale_color_manual(values = c("black", "darkgray")) +
  ylim(-.35, .35) +
  facet_wrap(~post_cw) +
  theme_bw() +
  theme(legend.position = "bottom",
        axis.text.x = element_text(color = "black", size = 10),
        axis.title.y = element_text(size = 12),
        legend.title = element_text(size = 12)) +
  labs(y = "Effect of Democratic Transition on Pr(Conflict)",
       x = "",
       color = "Years Since Transition (Pooled)",
       shape = "Years Since Transition (Pooled)")

# save plot
ggsave("figures/pooled-cold-war.pdf", width = 6.5, height = 5, units = "in", bg = "white")
